// display the convolved momentum distribution
// the input file format is mom str tot
// by Asia Sun, Jun. 23, 2022
// Last modified: Feb. 8, 2023

static int co = 0; // graph color
void primp(TGraph *g); // primp the graph
const string dir = "";
// C16-n c.s. on C12 target @239MeV/u
// const double csExp = 82.7639, dcsExp = 8.3707;

void momConv(){
  string file = dir+"mom.txt";
  string file0str = dir+"0d55p_strTConv.txt";
  string file0tot = dir+"0d55p_totConv.txt";
  string file1str = dir+"1s11p_strTConv.txt";
  string file1tot = dir+"1s11p_totConv.txt";
  TGraph *gtot = new TGraph(file.c_str(), "%lg %*lg %lg"); primp(gtot);
  TGraph *gstr = new TGraph(file.c_str(), "%lg %lg %*lg"); primp(gstr);
  TGraph *gtot0 = new TGraph(file0tot.c_str()); primp(gtot0);
  TGraph *gstr0 = new TGraph(file0str.c_str()); primp(gstr0);
  TGraph *gtot1 = new TGraph(file1tot.c_str()); primp(gtot1);
  TGraph *gstr1 = new TGraph(file1str.c_str()); primp(gstr1);
  TLegend *lg = new TLegend(0.6, 0.65, 0.89, 0.85); lg->SetBorderSize(0);
  // lg->AddEntry(gstr, "str");
  lg->AddEntry(gtot, "Total");
  // lg->AddEntry(gstr0, "str. 0d5/2@5+");
  lg->AddEntry(gtot0, "0d_{5/2}\\otimes5/2^{+}");
  // lg->AddEntry(gstr1, "str. 1s1/2@1+");
  lg->AddEntry(gtot1, "1s_{1/2}\\otimes1/2^{+}");

  // the exp. momentum distribution //
  TFile *f = new TFile("ppl.root");
  TH1F *hmomdisNET = (TH1F *)f->Get("hmomdisNET1");
  lg->AddEntry(hmomdisNET, "Exp. [38]");
  // // output the exp distribution to text //
  // const int nb = hmomdisNET->GetNbinsX();
  // ofstream fout("ppl.txt");
  // for(int i = 0; i < nb; i++){
  //   double x = hmomdisNET->GetBinCenter(i);
  //   double y = hmomdisNET->GetBinContent(i);
  //   double dy = hmomdisNET->GetBinError(i);
  //   fout << x << " " << y << " " << dy << endl;
  // } // end for over i
  // fout.close();


  // scaling //
  TGraph *gstrS = new TGraph(); primp(gstrS);
  const int n = gstr->GetN();
  double x, y;
  const double ss[] = {1.167*pow(16./15, 2), 0.734*pow(16./15, 2)};
  const double dm = gtot->GetMean() - hmomdisNET->GetMean() - 30.;
  // const double dm = 0.;
  cout << "gtot->GetMean(): " << gtot->GetMean() << endl;
  cout << "hmomdisNET->GetMean(): " << hmomdisNET->GetMean() << endl;
  cout << "dm: " << dm << endl;
  const double cc = gtot->Integral()/gstr->Integral();
  const double cce = hmomdisNET->Integral("width")/gtot->Integral(); // 320.; //
  // const double cce = 1.;
  cout << "gtot->Integral(): " << gtot->Integral() << endl;
  cout << "hmomdisNET->Integral(): " << hmomdisNET->Integral() << endl;
  cout << "cce: " << cce << endl;
  for(int i = 0; i < n; i++){
    gtot->GetPoint(i, x, y); gtot->SetPoint(i, x-dm, y*cce);
    gstr->GetPoint(i, x, y); gstrS->SetPoint(i, x-dm, y*cc*cce);
    gstr0->GetPoint(i, x, y); gstr0->SetPoint(i, x-dm, y*ss[0]*cce);
    gtot0->GetPoint(i, x, y); gtot0->SetPoint(i, x-dm, y*ss[0]*cce);
    gstr1->GetPoint(i, x, y); gstr1->SetPoint(i, x-dm, y*ss[1]*cce);
    gtot1->GetPoint(i, x, y); gtot1->SetPoint(i, x-dm, y*ss[1]*cce);
  } // end for over i
  // lg->AddEntry(gstrS, "str.-ScaledToTot");

  TCanvas *c = new TCanvas("c", "c", 800, 600);
  // gstrS->Draw();
  gtot->SetLineWidth(2); gtot->SetLineStyle(1);
  gtot->SetTitle("");
  TAxis *ax = gtot->GetXaxis(); TAxis *ay = gtot->GetYaxis();
  ax->SetTitleSize(0.05); ax->SetTitleOffset(0.91); ax->SetLabelSize(0.043);
  ay->SetTitleSize(0.05); ay->SetTitleOffset(0.98); ay->SetLabelSize(0.045);
  gtot->Draw();
  // gstr->Draw("same");
  gtot0->SetLineWidth(2); gtot0->SetLineStyle(2);
  gtot0->Draw("same");
  // gstr0->Draw("same");
  gtot1->SetLineWidth(2); gtot1->SetLineStyle(10);
  gtot1->Draw("same");
  // gstr1->Draw("same");
  gStyle->SetOptStat(0);
  hmomdisNET->Draw("same");
  gStyle->SetOptStat(0);

  lg->Draw();
  c->Print("/home/asia/Lab/cnok/cnok/docs/mom22/momNet.pdf");
} // end main function mom

// primp the graph
void primp(TGraph *g){
  TAxis *ax = g->GetXaxis(); ax->SetTitle("p_{c} [MeV/c]"); ax->CenterTitle();
  TAxis *ay = g->GetYaxis(); ay->SetTitle("d\\sigma_{-n} /dp_{c} [mb/(MeV/c)]");
  ay->CenterTitle(); ay->SetTitleOffset(1.30);
  co++;
  if(5 == co || 10 == co) co++;
  g->SetLineColor(co); g->SetMarkerColor(co);
} // end function primp
